This project aimed to look at the spatial variability of Symbiodinium clades C and D in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distribution of symbionts at scales ranging from location in the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch), colony color morph (brown vs. orange) and depth. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of stress-response potential in Kane’ohe Bay.
setwd("~/Symcap")
library(data.table)
library(devtools)
library(plyr)
library(reshape2)
library(popbio)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)
Coral_Data <- read.csv("Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path = "qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
target.ratios=c("C.D"),
fluor.norm = list(C=2.26827, D=0),
copy.number=list(C=33, D=3))
Mcap <- Mcap$result
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
Mcap <- Mcap[with(Mcap, order(Colony)), ]
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")
Symcap<-merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
Symcap$Mix <- factor(ifelse(Symcap$propC>Symcap$propD, ifelse(Symcap$propD!=0, "CD", "C"), ifelse(Symcap$propD>Symcap$propC, ifelse(Symcap$propC!=0, "DC", "D"), NA)), levels = c("C", "CD", "DC", "D"))
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
To account for the influence of tides, depth was adjusted according to the difference in sea level from the mean sea level on each collection day at 6-minute intervals. Mean sea level was obtained from NOAA tide tables for Moku o Lo’e.
JuneTide=read.csv("Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Station_1612480_tide_ht_20160701-20160731.csv")
Tide<-rbind(JuneTide, JulyTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)), format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2
Round6 <- function (times) {
x <- as.POSIXlt(times)
mins <- x$min
mult <- mins %/% 6
remain <- mins %% 6
if(remain > 3L) {
mult <- mult + 1
} else {
x$min <- 6 * mult
}
x <- as.POSIXct(x)
x
}
Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
merged<-merge(Symcap, Tide, by="Time.r", all.x=T)
merged$newDepth <- merged$Depth..m.- merged$TideHT
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Dom, Symcap$Reef.Area)
chisq.test(results)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 136.26, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##
## Slope Top
## C 0.7767857 0.3294574
## D 0.2232143 0.6705426
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Symbmiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
results=table(Symcap$Dom, Symcap$Color.Morph)
chisq.test(results)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 164.96, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##
## Brown Orange
## C 0.8896321 0.4103194
## D 0.1103679 0.5896806
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Color Morph", ylab = "Symbiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
results=table(Symcap$Mix, Symcap$Color.Morph)
chisq.test(results)
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 167.44, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##
## Brown Orange
## C 0.762541806 0.361179361
## CD 0.127090301 0.049140049
## DC 0.107023411 0.570024570
## D 0.003344482 0.019656020
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Color Morph", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Mix, Symcap$Reef.Area)
chisq.test(results)
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 138.97, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##
## Slope Top
## C 0.678571429 0.275193798
## CD 0.098214286 0.054263566
## DC 0.214285714 0.651162791
## D 0.008928571 0.019379845
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Reef Area", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)
results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 707, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##
## C D
## C 0.86635945 0.00000000
## CD 0.13364055 0.00000000
## DC 0.00000000 0.96703297
## D 0.00000000 0.03296703
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray 10", "gray 85", "gray 40", "gray100"), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray85", "gray40", "gray100"), inset = c(-.2, 0), xpd = NA)
propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
hist(propD, xlab = "Proportion of Clade D", ylab = "Number of Samples", main = "", col = "gray75")
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Color.Morph, Symcap$Reef.Area)
chisq.test(results)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 81.109, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##
## Slope Top
## Brown 0.5523385 0.2015504
## Orange 0.4476615 0.7984496
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Color Morph Proportion")
legend("topright", legend=c("Brown", "Orange"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
results=table(Symcap$Color.Morph, Symcap$Reef.ID)
chisq.test(results)
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 57.007, df = 25, p-value = 0.0002667
prop.table(results, margin = 2)
##
## 10 13 14 18 19 20
## Brown 0.4000000 0.4000000 0.5333333 0.2333333 0.3666667 0.6285714
## Orange 0.6000000 0.6000000 0.4666667 0.7666667 0.6333333 0.3714286
##
## 21 25 26 30 37 38
## Brown 0.4333333 0.6000000 0.4242424 0.1333333 0.8000000 0.2000000
## Orange 0.5666667 0.4000000 0.5757576 0.8666667 0.2000000 0.8000000
##
## 42 46 5 Deep F1-46 F2-25
## Brown 0.5142857 0.5000000 0.3666667 0.5172414 0.4000000 0.6000000
## Orange 0.4857143 0.5000000 0.6333333 0.4827586 0.6000000 0.4000000
##
## F3-18 F4-34 F5-34 F6-Lilipuna F7-18 F8-10
## Brown 0.2500000 0.4400000 0.2500000 0.5000000 0.2000000 0.4000000
## Orange 0.7500000 0.5600000 0.7500000 0.5000000 0.8000000 0.6000000
##
## F9-5 HIMB
## Brown 0.4000000 0.6000000
## Orange 0.6000000 0.4000000
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
Dom1 <- subset(merged, !is.na(newDepth) & !is.na(Dominant))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
logi.hist.plot(Dom1$newDepth, Dom1$Dominant, boxp = FALSE, type = "hist", col="gray", xlabel = "Depth (m)", ylabel = "", ylabel2 = "")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "C D", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of clade C Symbiont", line = 3, cex = 1)
merged$DepthInt <- cut(merged$Depth..m., breaks = 0:13)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.2, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Dominant~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Dominant Symbiont Proportion")
Bars indicate relative frequency of clade C vs. D dominance at 1m depth intervals when pooling all samples collected.
merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
merged$Mixture2 <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 0, 1)
results=glm(Mixture~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Mixture
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 589 817.24
## newDepth 1 75.664 588 741.57 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
results=table(merged$Mixture2, merged$DepthInt)
results
##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 172 43 30 34 18 11 3 2 5 2 1
## 1 66 41 60 99 46 20 10 17 12 6 6
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Mix", "Non-Mix"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.23, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Mixture~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Mixture Proportion")
merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
results=glm(Color~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
Color <- subset(merged, !is.na(newDepth) & !is.na(Color))
logi.hist.plot(independ = Color$newDepth, depend = Color$Color, type = "hist", boxp = FALSE, ylabel = "", col="gray", ylabel2 = "", xlabel = "Depth (m)")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "Brown Orange", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Color Morph Proportion")
Bars indicate relative frequency of Brown vs. Orange color morph dominance at 1m depth intervals when pooling all samples collected.
merged$Reef.Area <- ifelse(merged$Reef.Area!="Top", yes = "Slope", no = "Top")
table(merged$Dom, merged$Color.Morph, merged$Reef.Area)
## , , = Slope
##
##
## Brown Orange
## C 226 122
## D 21 79
##
## , , = Top
##
##
## Brown Orange
## C 40 45
## D 12 161
model=aov(Dominant~Color.Morph*Reef.Area, data = merged)
Anova(model, type = 2)
## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## Color.Morph 21.329 1 134.208 < 2.2e-16 ***
## Reef.Area 14.489 1 91.168 < 2.2e-16 ***
## Color.Morph:Reef.Area 1.780 1 11.201 0.0008612 ***
## Residuals 111.565 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Because and interactive effect of reef area and color morph was observed and slope is indicative of a depth gradient, the interaction between depth and color morph was tested.
model1=lm(Dominant~Color.Morph*newDepth, data = merged)
Anova(model1, type = 2)
## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## Color.Morph 24.863 1 154.391 < 2.2e-16 ***
## newDepth 12.236 1 75.981 < 2.2e-16 ***
## Color.Morph:newDepth 2.592 1 16.097 6.798e-05 ***
## Residuals 94.369 586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model3=aov(Dominant2~Depth..m.*Reef.Type, data = merged)
Anova(model3, type = 2)
## Anova Table (Type II tests)
##
## Response: Dominant2
## Sum Sq Df F value Pr(>F)
## Depth..m. 24.214 1 120.3238 < 2.2e-16 ***
## Reef.Type 0.052 1 0.2569 0.612431
## Depth..m.:Reef.Type 1.727 1 8.5815 0.003506 **
## Residuals 141.272 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Dominant2~Depth..m., family = "binomial", data = df)
newdata <- list(Depth..m.=seq(0,12,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,12,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Dominant2~Depth..m., family = "binomial", data = df)
newdata <- list(Depth..m.=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,12,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade C Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)
While brown was always C-dominated, orange was observed to switch from D to C dominance. The depth threshold at which this switch occurs was analyzed here.
threshdepth <- function(color) {
df <- subset(merged, Color.Morph==color)
plot(df$Dominant2~df$newDepth, xlab="Depth (m)", ylab = "Probability of Clade C Symbiont",
main=color)
abline(h = 0.5, lty=2)
results=glm(Dominant2~newDepth, family = "binomial", data = df)
pval <- data.frame(coef(summary(results)))$`Pr...z..`[2]
mtext(side=3, text=pval)
newdata <- list(newDepth=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted ~ seq(0,12,0.01))
thresh <- ifelse(pval < 0.05,
newdata$newDepth[which(diff(sign(fitted - 0.5))!=0)], NA)
return(thresh)
}
sapply(levels(merged$Color.Morph), FUN=threshdepth)
## Brown Orange
## NA 2.7
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,12,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,12,0.01), col="sienna", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade C Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.22, 0), xpd = NA)
model2=aov(Color~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)
## Anova Table (Type II tests)
##
## Response: Color
## Sum Sq Df F value Pr(>F)
## newDepth 6.180 1 27.6214 2.069e-07 ***
## Reef.Type 0.132 1 0.5908 0.442436
## newDepth:Reef.Type 1.651 1 7.3813 0.006785 **
## Residuals 131.109 586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,12,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,12,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)
KB <- c(21.46087401, -157.809907)
KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, col=153, pch=21, bg="#FF7F50", lwd=2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
par(new=T, mar=c(11.8,19,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
par(mfrow=c(3,1))
merged$DepthInt <- cut(merged$Depth..m., breaks = 0:13)
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(2, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(2.1, 4, 2, 6))
results=glm(Dominant~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="Probability of D-Dominance")
merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(3, 4, 1, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "Probability of Orange-Dominance",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(3.1, 4, 1, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="")
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 0, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="Depth (m)", ylab="Probabilty of D-Dominance", axisnames=FALSE, xaxs = "i", yaxs = "i")
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0, 11, 0.01), col="sienna", lwd=3)
#legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.13, 0), xpd = NA)
results=table(merged$Dominant2, merged$Color.Morph)
chisq.test(results)
prop.table(results, margin = 2)
par(new=T, mar=c(10, 10, .5, 6.3))
barplot(prop.table(results, margin = 2), col = c(alpha("red", 0.35), alpha("blue", 0.35)), xlab = "", ylab = "", yaxt = 'n')
#legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.35), alpha("red", 0.35)), inset = c(0, 0), xpd = NA)